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Abstract. We have developed finite difference codes based on the Yin- Yang grid for the 
geodynamo simulation and the mantle convection simulation. The Yin- Yang grid is a kind of 
spherical overset grid that is composed of two identical component grids. The intrinsic simplicity 
of the mesh configuration of the Yin- Yang grid enables us to develop highly optimized simulation 
codes on massively parallel supercomputers. The Yin- Yang geodynamo code has achieved 15.2 
Tflops with 4096 processors on the Earth Simulator. This represents 46% of the theoretical 
peak performance. The Yin- Yang mantle code has enabled us to carry out mantle convection 
simulations in realistic regimes with a Rayleigh number of 10 7 including strongly temperature- 
dependent viscosity with spatial contrast up to 10 6 . 



1. Introduction 

The Earth (radius r = 6400km) is composed of three spherical layers; the inner core (r = 
1200km), the outer core (r = 3500km), and the mantle. Computer simulations of the Earth's 
interior need efficient spatial discretization methods in the spherical shell geometry. To achieve 
high sustained performance on massively parallel supercomputer such as the Earth Simulator, 
spatially localized discretization methods rather than spectral methods are desirable. Recently, 
we proposed a new spherical grid system, the "Yin- Yang grid," for geophysical simulations. 
Because there is no grid mesh that is orthogonal over the entire spherical surface and, at the 
same time, free of coordinate singularity or grid convergence, we have chosen an overset grid 
approach. A spherical surface is decomposed into two identical subregions. The decomposition 
(or dissection) enables us to cover each subregion by a grid system that is individually orthogonal 
and singularity-free. Each component grid in this Yin- Yang grid is a low latitude component of 
the usual latitude-longitude grid on the spherical polar coordinates (90 degree about the equator 
and 270 degree in the longitude). Therefore, the grid spacing is quasi- uniform and the metric 
tensors are simple and analytically known. One can directly apply mathematical and numerical 
resources that have been written in the spherical polar coordinates or latitude-longitude grid. 
Since the two component grids are identical and they are combined in a complementary way, 
various routines of the code can be recycled twice for each component grid at every simulation 
time step. We have developed finite difference codes based on the Yin- Yang grid for (i) the 
geodynamo simulation in the outer core, and (ii) the mantle convection simulation. 

In general, a dissection of a computational domain generates internal borders or internal 
boundaries between the subregions. In the overset grid methodology PP, the subregions are 



permitted to partially overlap one another on their borders. The overset grid is also called as 
overlaid grid, or composite overlapping grid, or Chimera grid |2j. The validity and importance 
of the overset approach in the aerodynamical calculations was pointed out by Steger [Hj. Since 
then this method is widely used in this field. It is now one of the most important grid techniques 
in the computational aerodynamics. 

In the computational geosciences, the idea of the overset grid approach appeared rather early. 
Phillips proposed a kind of composite grid in 1950's to solve partial differential equations on 
a hemisphere, in which the high latitude region of the latitude-longitude grid is "capped" by 
another grid system that is constructed by a stereographic projection to a plane on the north 
pole |IJ|!)1 After a long intermission, the overset grid method seems to attract growing interest 
in geoscience these days. The "cubed sphere" 7 is an overset grid that covers a spherical surface 
with six component grids that correspond to six faces of a cube. The "cubed sphere" is recently 
applied to the mantle convection simulation 8 . In the atmospheric research, other kind of 
spherical overset grid is used in a global circulation model |5], in which the spherical surface 
is covered by two component grids — improved stereographic projection grids — in northern and 
souther hemispheres that overlap in the equator. 

Among indefinite variations of spherical overset grid systems, what is the simplest one? In 
general, the structure of a spherical overset grid is largely determined by the number of divided 
pieces of the sphere n (> 2). Here we consider the minimum case of n = 2, i.e., the spherical 
dissections by two pieces. One can divide a sphere into two parts, for example, by cutting along 
a small circle at any latitude. We concentrate on a special class of n = 2 dissections in which 
the two pieces are geometrically identical, i.e., they have exactly same size and shape. Another 
condition we impose here to maximize the simplicity is the symmetry of the piece. It should 
have two fold symmetry in two perpendicular directions; up-down and right-left. Here we call 
this special class of dissections as yin-yang dissection of a sphere. 

A trivial example of the yin-yang dissection is obtained by cutting along the equator or any 
great circle, producing two hemispheres. 
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Figure 1. An example of yin-yang dissection 
of a sphere: A sphere is divided into two 
identical pieces, with same shape and size. 
Each piece has two fold symmetry; up- 
down and right-left. They are combined 
in a complemental way to cover a spherical 
surface. The two identical pieces of the yin- 
yang dissection is transformed each other by 
two successive rotations, or one rotation. 



Other yin-yang dissections are obtained by modifying the cut curve from the great circle. Let 
Syin be a piece of a sphere S with radius r = y/2. We should keep the surface area of S y i n being 
2-Kr 2 , just a half of S's surface. An example of S y i n is shown in the upper left panel in Fig. ^ 



The border curve of S y i n passes through the following four points on the sphere; point A at 
(x,y,z) = (0,-1, +1), B at (0,+l,+l), C at (0,+l,-l) and D at (0,-1,-1). The curve AB, 
between A and B, is arbitrarily as long as it is symmetric about the y = plane. Other three 
curves, BC, CD, and DA, are uniquely constructed from the curve AB as follows: The curve 
BC is a copy of AB followed by two successive rotations, first 180 degree about the z axis, then 
90 degree about the x axis. The curve CD is the mirror image of AB about z = plane. The 
curve DA is the mirror image of BC about y = plane. Prom this definition of the border curve 
ABCD, it is obvious that the surface area of S y i n is just a half of that of the sphere S. Now 
we make a copy of S y i n and call it S yang which is rotated for 180° around z-axis. (See lower left 
panel of Fig. ^) Then, rotate it again, but this time for 90° degree around shown in 

the lower right panel. Then the original piece S y i n (the upper left) and the rotated copy S yany 
(the lower right) can be combined, and they just cover the original sphere S as shown in the 
upper right in this figure. This is an constructive illustration of the yin-yang dissection of a 
sphere. 

Since the initial curve AB was arbitrarily, it is obvious that there are indefinite variations of 
the yin-yang dissection of the sphere S. 



2. Yin- Yang grids 




Figure 2. A dissection of a sphere into 
two identical pieces — Yin and Yang — with a 
partial overlap. The thick and thin curves 
are the borders of Yin and Yang piece, 
respectively, the thick curve (Yin's border) is 
always located in either constant latitudes or 
constant longitudes. The Yin (Yang) piece is 
a rectangle in the computational (9, (j)) space 
of the Yin (Yang) grid. 



The overset grid methodology gives us a freedom to design the shape of the component grid as 
long as the grids has minimum overlap one another pQ. Therefore, we can take the component 
grid as a rectangle in the computational (0,4>) space. Fig. El shows a spherical dissection by 
two identical pieces with partial overlap. The Yin piece is surrounded by a thick red curve and 
Yang piece is surrounded by a thin curve. Note that the northern and southern borders of the 
Yin piece are located in constant latitudes and the western and eastern borders are located in 
constant longitudes. In other words, the Yin piece in Fig. [21 is a rectangle in the (6,(j)) space 
of the Yin's spherical coordinates, and therefore, the Yang piece is also (the same) rectangle 



in Yang's coordinates that is perpendicular to the Yin's. The Yin- Yang grid based on this 
partially overlapped spherical dissection is shown in Fig. El Here, each component grid spans 
the subregion S y defined by 

S y := {9, (/)}, I 9 - tt/2 I < tt/4 + 5, | 4> | < 3^/4 + 5, (1) 

with a small buffer 5 which is necessary to keep the minimum overlap between Yin and Yang. 
Note that in the simulation code, one subroutine for the fluid solver, for instance, can be recycled 
twice because the grid distribution is exactly the same for the Yin and Yang. 




Figure 3. A Yin- Yang grid based on the yin-yang dissection with partial overlap shown in 
Fig. |21 Each component grid is rectangle in the computational (#,</>) space. 

The Yin and Yang are converted each other by a rotation. The Yin's cartesian coordinates 
x- 1 for i = 1,2,3 and that Yang's coordinates x\ are related by 

xl = M ijX ] for i,j = 1,2,3 (2) 

where M\\ = — 1, M23 = M32 = 1, and Mu = for other components. Note that the matrix M 
satisfies 

M = M l = M'\ (3) 

which indicates a complemental relation between the Yin and Yang. The coordinate 
transformation from Yin to Yang is mathematically the same as that from Yang to Yin. This 
enables us to make only one, instead of two, subroutines that involve any data transformation 
between Yin and Yang, which is required in the mutual interpolation for the internal boundary 
condition on the overset grid borders. 

The transformation formula of any vector components v = (v r , vq^v^) between Yin and Yang 
is given by 
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with the transformation matrix 
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The inverse transformation is given by the same matrix; P^ 1 = P, which is another reflection 
of the complemental nature between the Yin and Yang. 

Another merit of the Yin- Yang grid resides in the fact that the component grid is nothing 
but (a part of ) the latitude-longitude grid. We can directly deal with the equations to be solved 
with the usual spherical polar coordinates. The analytical form of metric tensors are familiar in 
the spherical coordinates. We can directly code the basic equations in the program as they are 
formulated in the spherical coordinates. We can make use of various resources of mathematical 
formulas, program libraries, and other tools that have been developed in the spherical polar 
coordinates. 

In order to illustrate the programing strategy in the Yin- Yang method, let us consider a 
two-dimensional fluid problem on a sphere S. Suppose that two components of the flow velocity 
v = (vQ,V(h) and the pressure p are written in vel_t, vel_p, and press in a Fortran 90/95 
program. They can be combined into one structure or "type" in Fortran 90/95 as 

type fluid. 

real(DP), dimension (NT, NP) :: vel_t, vel_p, press 
end type fluid_ 

where NT, NP are the grid size integers in 6 and directions in the subregion S y of eq. Using 
this structured type, we declare two variables for the fluid; one is for Yin and another for Yang: 

type(fluid_) : : fluid_yin, fluid_yang 

Then, we call a fluid solver subroutine, here named navier_stokes_solver, that numerically 
solves the Navier-Stokes equation in the spherical coordinates in the subregion S y : 

call navier_stokes_solver (f luid_yin) 
call navier_stokes_solver (f luid_yang) 

The first call of navier_stokes_solver solves the fluid motion in the S y region defined in the 
Yin's spherical coordinates and the second call is for the same region S y defined in the Yang's 
coordinates. But in the program code, we do not have to distinguish the two S y regions since 
the basic equations, numerical grid distribution, and therefore, all numerical tasks are identical 
in the computational space. For a rotating fluid problem with a constant angular velocity ft, we 
have the Coriolis force term in the Navier-Stokes equation that seems to break the symmetry 
between the Yin grid and Yang grid, but it is still possible to write the equation in exactly the 
same form for the Yin and Yang grids by explicitly writing three components of angular velocity 
in the Coriolis force term 2v x CI in the subroutine. Then, we call the routine with the angular 
velocity vector in each grid (Yin or Yang) as the second argument: 

call navier_stokes_solver (f luid_yin, omega_yin) 
call navier_stokes_solver (f luid_yang, omega_yang) 

where omega_yin and omega_yang are again structured variables that hold three components of 
the fi vector: For example, omega_yin holds three components of cartesian vector components 
in the Yin grid (O", 0™, 0") = (0, 0, Q), and omega_yang holds (fl e x , Q e y , J2«) = (0, 0). 

Our experience tells that it is easy to convert an existing latitude-longitude based program 
into a Yin- Yang based program since there are many shared routines between them. In addition 
to that the size of the code as well as its complexity is drastically reduced by the code conversion 
because we can remove routines that are designed to resolve the pole problems on the latitude- 
longitude grid. 
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3. Application to the mantle convection simulation 

3.1. Simulation model 

We applied the Yin- Yang grid described in the previous secion for the mantle convection 
simulation. The details of the adopted numerical methods and benchmark tests can be found 
in |ng. 

We model the mantle convection as a thermal convection of a Boussinesq fluid with 
infinite Prandtl number heated from bottom of a spherical shell The ratio of the inner 

radius (r = tq) and the outer radius (r = n) is 0.55. The normalization factors for the 
non-dimensionalization of the length, velocity, time and temperature are d, k/d, d 2 /k and 
AT 

— Tbot ~ T( p, respectively, where d is the thickness of the shell, k the thermal diffusivity, 
and Tfcot and Tt op are the temperatures on the bottom and top surfaces. The hat stands 
for dimensional quantity. The non-dimensional equations of mass, momentum, and energy 
conservation governing the thermal convection are, 

V • v = 0, (6) 

= -Vp + V • (jib) + RaTr, (7) 
dT 

- = V 2 T-vVT + ff, (8) 

where v is the velocity vector, p pressure, ji viscosity, T temperature, t time, e strain-rate tensor, 
and f is the unit vector in the r-direction. The Rayleigh number is defined by 

where p is the density, g the gravitational acceleration, and a is the thermal expansivity. Most 
of the heat for Earth's mantle comes from a combination of radioactive decay of isotopes and 
secular cooling of the mantle. The internal heating is defined by 

KCpAl 

where Q is the internal heating rate per unit mass, and c p is the specific heat at constant 
pressure. 

According to the laboratory experiments on silicate rock deformation, the viscosity of the 
Earth's mantle depends on various parameters such as temperature, pressure, stress, and so 
on |12| . Among them, temperature dependence is the most dominant factor. Here we assume 
that viscosity \i depends only on temperature; 

/i(T) = exp[-£(T-T 6ot )]. (11) 

The viscosity contrast across the spherical shell is defined by 7 M = fi(T top ) / /i(?6 ot ) = exp(E). 
The mechanical boundary conditions at the top and bottom surface are immpermiable and 
stress-free. The boundary conditions for T are fixed; Tf, t = 1 and Tt op = 0. 



3.2. Steady state convection 

The thermal convection in the spherical shell with infinite Prandtl number has two stable 
solutions with polyhedral symmetry when the Rayleigh number is low. The two solutions 
are found by linear theory and confirmed by numerical simulations One solution is a 

convection with the tetrahedral symmetry which has four upwellings; the other has the cubic 




Figure 4. The iso-surface of the residual temperature 5T (the deviation from horizontally 
averaged temperature at each depth) started from the initial condition of (a) the tetrahedral 
and (b) the cubic symmetries. The Rayleigh number is Ra = 10 4 . Blue and Yellow iso-surfaces 
indicate 5T = —0.125 and 5T = 0.150, respectively. Red spheres indicate the bottom of the 
mantle with fixed temperature. 



symmetry with six upwellings. To confirm these symmetric solutions and their stabilities, we 
performed two simulations with different initial conditions of temperature field; T(r, 9, (ft) = 
Tcond(r) +T prtb (r, 9, (ft), where T cond (r) is the purely conductive profile, V 2 T cond (r) = 0, with the 
thermal boundary conditions given above. The perturbation term T pr tb(r,9,(ft) is given by, 

T prtb (r,9,(ft) =0.1 Y 3 2 (9,(ft) sin Tr(r-r ), (12) 

for the tetrahedral symmetric solution, and 

T prtb (r, 9, (ft) = 0.1 {y 4 °(0, (ft) + ^Y A \9, 0) j sinvr(r - r„), (13) 

for the cubic symmetric solution, where Y( m (9, (ft) is the normalized spherical harmonic functions 
of degree I and order ra. Fig. shows the steady state convection pattern with the tetrahedral 
and cubic symmetries. We have performed benchmark tests with previously reported numerical 
mantle convection codes that employed various numerical schemes. In spite of the differences of 
the discretization methods, numerical techniques, and number of grid points among the codes, 
we found that the calculated values such as the Nusselt number obtained by our Yin- Yang 
mantle code agree well with previous calculations within a few percent. 

3.3. Time- dependent convection 

The Earth's mantle is obviously in a time-dependent convection under a very high Rayleigh 
number (Ra > 10 6 ) and with internal heating (H < 20). When Ra = 10 5 , the convection pattern 
becomes weakly time-dependent, and the geometrical symmetry is broken. Fig. 03 shows the 
thermal structures of the mantle convection when Ra = 10 which is characteristic of the Earth's 
mantle. Without internal heating, the thermal structure is strongly time-dependent, driven by 
narrow, cylindrical upwelling (hot) plumes surrounding by a network of long downwelling (cold) 
sheets (Fig. EH)- This feature is in contrast with the convective feature at low Rayleigh number 
(Ra < 10 5 ) where the convection is nearly steady state (Fig. QJ. On the other hand, when the 




Figure 5. The iso-surface of the temperature T and the residual temperature 5T for the cases 
of (a) H = and (b) H = 20. The Rayleigh number is Ra = 10 7 . Iso-surfaces on the half 
spherical shell indicate the temperature (see color bars). Blue and Yellow iso-surfaces indicate 
5T = ±0.1, respectively. Red spheres indicate the bottom of the mantle with fixed temperature. 



internal heating is taken into account {H = 20), the convective feature is dominated by the 
short-wavelength structure with numerous quasi-cylindrical downwellings spaced relatively close 
together. The downwellings are surrounded by a broad and diffuse upwelling of hotter fluid 
(Fig. 0a). We have found that internal heating has a strong influence on the scale and structure 
of the mantle convection, especially on the shape of downwellings. 

The convection pattern is also drastically changed by taking the viscosity variation into 
account. Fig. shows the thermal structures of the mantle convection with temperature- 
dependent viscosity at Ra = 10 7 and H = 0. When the temperature dependence of viscosity is 
rather moderate (the viscosity contrast across the convecting shell 7^ is 10 3 -10 4 ), the convection 




Figure 6. The iso-surface of the residual temperature for the cases of (a) ^y mu = 10 4 (E = 9.210) 
and (b) 7^ = 10 6 [E = 13.816). Blue and Yellow iso-surfaces indicate (a) 6T = ±0.25 and 
ST = ±0.10, respectively. Red spheres indicate the bottom of the mantle with fixed temperature. 



has long-wavelength thermal structure with a mobile, stiff layer, or, "sluggish-lid" along the cold 
top surface of the mantle. When Ra = 10 7 and 7^ = 10 4 , the convection pattern comes to be 
dominated by the degree-one pattern; the one cell structure that consists of a pair of cylindrical 
downwelling plume and cylindrical upwelling plume (Fig. EH)- On the other hand, the convective 
flow pattern that belongs to the "stagnant-lid" regime emerges when 7^ > 10 5 . The stagnant-lid, 
which is an immobile, stiff layer, prevents the heat flux through the top boundary and leads to 
a small temperature difference in the mantle below the lid. Convection under the stagnant-lid is 
characterized by numerous, small-scale cylindrical plumes surroundings sheet-like downwelling 
(Fig. m>). We have found that the variable viscosity with temperature dependence induces 
drastic effects on the mantle convection pattern. 



4. Application to geodynamo simulation 

The magnetic compass points to the north since the Earth is surrounded by a dipolar magnetic 
field. It is broadly accepted that the geomagnetic field is generated by a self-excited electric 
current in the Earth's core, The inner core is iron in solid state, and the outer core is also iron 
but in liquid state due to the high temperature of the planetary interior. The electrical current 
is generated by magnetohydrodynamic (MHD) dynamo action — the energy conversion process 
from flow energy into magnetic energy — of the liquid iron in the outer core. In the last decade, 
computer simulation has emerged as a central research method for geodynamo study |14j . 

In this section, we show the application of the Yin- Yang grid to the geodynamo simulation 
with a special emphasize on the code parallelization and sustained performance achieved by the 
Earth Simulator. We consider a spherical shell vessel bounded by two concentric spheres. The 
inner sphere of radius r = r j denotes the inner core and the outer sphere of r = r Q denotes the 
core-mantle boundary. An electrically conducting fluid is confined in this shell region. Both 
the inner and outer spherical boundaries rotate with a constant angular velocity CI. We use a 
rotating frame of reference with the same angular velocity. There is a central gravity force in 
the direction of the center of the spheres. The temperatures of both the inner and outer spheres 
are fixed; hot (inner) and cold (outer). When the temperature difference is sufficiently large, a 
convection motion starts when a random temperature perturbation is imposed at the beginning 
of the calculation. At the same time an infinitesimally small, random "seed" of the magnetic 
field is given. 

The system is described by the following normalized MHD equations: 

|-V.f, (14) 
df 1 

— = -V- (vf) - Vp + j x B + / og + 2pv x fi + ^(V 2 v + -V(V • v)), (15) 

§ = -v • Vp - 7 pV • v + (7 - l)KV 2 T + (7 - l)r?j 2 + (7 - 1)$, (16) 
at 

^ = vxB + ,V 2 A, (17) 

with 



p = pT, B = V x A, j = V x B, g = -g /r 2 r, 

V-B = 0, $ = 2 M (e^-±(V-v) 2 ), e ^. = i(f| + §|). (18) 

Here the mass density p, pressure p, mass flux density f , magnetic field's vector potential A are 
the basic variables in the simulation. Other quantities; magnetic field B, electric current density 
j, and electric field E are treated as subsidiary fields. The ratio of the specific heat 7, viscosity 



Table 1. Specifications of the Earth Simulator. 



Peak performance of arithmetic processor (AP) 8 Gflops 

Number of AP in a processor node (PN) 8 

Total number of PN 640 

Total number of AP 8 AP x 640 PN = 5120 

Shared memory size of PN 16 GB 

Total peak performance 8 Gflops x 5120 AP = 40Tflops 

Total main memory 10 TB 

Inter-node data transfer rate 12.3 GB/s x 2 



fj,, thermal conductivity K and electrical resistivity rj are assumed to be constant. The vector 
g is the gravity acceleration and f is the radial unit vector; g$ is a constant. We normalize the 
quantities as follows: The radius of the outer sphere r Q = 1; the temperature of the outer sphere 
T(l) = 1; and the mass density at the outer sphere p{l) = 1. The temperature on the inner and 
outer spheres are fixed. The boundary condition for the velocity is rigid; 

v = 0, atr = rj,l. (19) 

The boundary condition for the magnetic field is given by 

B e = B (f) = 0, atr = r;,l. (20) 

We will consider the improvement of this rather artificial boundary condition into more realistic 
one in the end of this section. The spatial derivatives in the above equations are discretized by 
the second-order central finite difference method on the Yin- Yang grid. The fourth-order Runge- 
Kutta method is used for the temporal integration. Initially, both the convection energy and the 
magnetic energy are negligibly small. For geodynamo study, it is necessary to follow the time 
development of the MHD system until the thermal convection flow and the dynamo-generated 
magnetic field are both sufficiently developed and saturated. 

We developed this Yin- Yang based geodynamo simulation code for the Earth Simulator by 
converting our previous geodynamo code, which was based on the traditional latitude-longitude 
grid, into the Yin- Yang grid. We have found that the code conversion from our previous latitude- 
longitude based code into the new Yin- Yang based code is straightforward and rather easy. Our 
experience with the rapid and easy conversion from latitude-longitude code into Yin- Yang code 
would be encouraging for others who have already developed codes that are based on latitude- 
longitude grids in the spherical coordinates, and who are bothered by numerical problems and 
inefficiency caused by the pole singularity. We would like to suggest that they try the Yin- Yang 
grid. 

Since the Yin grid and Yang grid are identical, dividing the whole computational domain into 
a Yin grid part and a Yang grid part is not only natural but also efficient for parallel processing. 
In addition to this Yin-and-Yang division, further domain decomposition within each grid is 
applied to for the massively parallel computation on the Earth Simulator. 

The Earth Simulator, whose hardware specifications are summarized in Table ^ has three 
different levels of parallelization: Vector processing in each arithmetic processor (AP); shared- 
memory parallelization by 8 APs in each processor node (PN); and distributed-memory 
parallelization by PNs. 



In our Yin- Yang dynamo code, we apply vectorization in the radial dimension of the three- 
dimensional (3D) arrays for physical variables. The radial grid size is 255 or 511, which is just 
below the size (or doubled size) of the vector register of the Earth Simulator (256) to avoid bank 
conflicts in the memory. We use MPI both for the inter-node (distributed memory) parallel 
processing and for the intra-node (shared memory) parallel processing. This approach is called 
"flat-MPI" parallelization. 

As we mentioned above, we first divide the whole computational domain into two identical 
parts that correspond to the Yin grid and Yang grid shown in Fig. |3fa). (Therefore, the total 
number of processes is always even.) For further parallelization within each component grid, we 
applied the two-dimensional decomposition in the horizontal space, colatitude 9 and longitude 
(j>. More details on the parallelization of this code is described in |15j . 

The best performance of the Yin- Yang geodynamo code with the flat MPI parallelization is 
15.2 Tflops. This performance is achieved by 4096 processors (512 nodes) with the total grid size 
of 511(radial) x 5 14 (latitudinal) x 1538(longitudinal) x 2(Yin and Yang). Since the theoretical 
peak performance of 4096 processors is 4096 x 8 Gflops = 32.8Tflops, we have achieved 46% 
of peak performance in this case. The average vector length is 251.6, and the vector operation 
ratio is 99%. The high performance of the Yin- Yang dynamo code is a direct consequence of 
the simple and symmetric configuration design of the Yin- Yang grid: It makes it possible to 
minimize the communication time (10%) between the processes in the horizontal directions, and 
enables optimum vector processing (with 99% of operation ratio) in the radial direction in each 
process. 

Before concluding this section, we briefly describe our recent improvement of the Yin- Yang 
geodynamo code. We have improved the boundary condition denoted by eq. (|2T?|) of the magnetic 
field into more realistic one, i.e., so called vacuum boundary condition. In this boundary 
condition, the magnetic field generated by the MHD dynamo in the outer core (r < 1) is 
smoothly connected to the magnetic field of the outer region r > 1 that is assumed to be an 
insulator; 

VxB„ = 0, for r > 1. (21) 
Therefore, the B„ is written by a scalar function ip, 
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-Vip, for r > 1, 



where, from V • Bt, = 0, ip satisfies the potential equation 
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ip = 0, for r > 1 



(22) 



(23) 



The boundary condition of ip at r = 1 is given by — Vip(r = 1) = B r {r = 1), where B r (r = 1) 
is determined from the dynamo region (r < 1). Other component of the magnetic field at the 
surface Bg(r = 1) and B^r = 1) are determined by the solution of eq. (|2*3*|) . In order to solve 
this boundary value problem, we first apply a coordinate transformation of r. 



r — > £ = 1/r. 

The equation (|23j) is converted into the following form 
d 2 13/ J\ Id 2 
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ip = 0, for < C < 1 



(24) 



(25) 



The problem to solve eq. (|23|) outside a unit sphere r > 1 is now converted into the problem to 
solve eq. (|25[) inside a unit sphere C < 1. The boundary condition of ip at the origin £ = is 
given by ^(C = 0) = since ip(r = oo) = 0. 




Figure 7. The Yin- Yang multigrid method for the solution of the vacuum magnetic field 
potential tp. The full approximation storage algorithm is used. The horizontal boundary values 
of the Yin and Yang grids for the overset are determined by the mutual interpolation (white 
arrows) at the every grid level in the V-cycle of the multigrid method (gray arrows). 



To solve eq. (j25j) , we apply the multigrid method |16j , which is practically the optimal way to 
solve this kind of boundary value problem. The base grid system is the Yin- Yang grid defined 
in the full spherical region including the origin. See Fig. [7| We adopt the full approximation 
storage algorithm of the multigird method. The Jacobi method is used as the smoother. The 
V-cycle is repeated for a couple of times until we get the convergence. The internal boundary 
condition of each component grid (Yin and Yang) are set by mutual bi-cubic interpolation at 
every grid level as indicated by white arrows in Fig. [Tj Although, the code is not parallelized 
yet, its flat-MPI parallelization will be straightforward. We have combined this non-parallelized 
Yin- Yang multigrid solver of the vacuum potential tjj with the non-parallel version of the Yin- 
Yang geodynamo code. We have found that the vacuum field condition has been successfully 
implemented by this multigrid potential solver with almost the same computational cost (CPU 
time) as with the MHD solver part. This is a very promising result for further development. 

5. Summary 

We have developed a new spherical overset grid, "Yin- Yang grid", for geophysical simulations. 
The Yin- Yang grid is constructed from a dissection of a sphere into two identical and 
complemental pieces. Among various possible overset grids over a sphere, we believe that the 
Yin- Yang grid is the simplest and the most powerful especially on massively parallel computers 
from the following reasons: 

• It is an orthogonal system, since it is a part of the latitude-longitude grid. 

• The grid spacing is quasi-uniform, since we picked up only the low latitude region of the 
latitude- longitude grid. 



• The metric tensors are simple and analytically known, since it is denned based on the 
spherical polar coordinates. 

• Routines for the fluid (or MHD) solver can be recycled twice, since Yin and Yang are 
identical. 

• Routines for mutual interpolations of the overset grid borders can also be recycled twice, 
since Yin and Yang are complemental. 

• Parallelization is easy and efficient, since the domain decomposition is straightforward. 

We have developed finite difference codes of the geodynamo simulation and the mantle 
convection simulation on the Yin- Yang grid. The Yin- Yang geodynamo code has achieved 15.2 
Tflops with 4096 processors on the Earth Simulator. This represents 46% of the theoretical 
peak performance. By the Yin- Yang mantle code, we can carry out realistic mantle convection 
simulations under the Rayleigh number of 10 7 , including strongly temperature-dependent 
viscosity whose contrast reaches upto 10 6 . 

In the Earth Simulator Center, the Yin- Yang grid is also applied to advanced general 
circulation modes of the atmosphere and ocean [T71 HH1 IT5| . 



Acknowledgments 

The authors would like to thank Prof. Tetsuya Sato, the director-genenal of the Earth Simulator 
Center, for instructive discussions and Dr. Masanori Kameyama for useful comments on the 
application of the multigrid method to the Yin- Yang grid. All the simulations were performed 
on the Earth Simulator. 



References 

[1] Chesshire G and Henshaw W D 1990 Composite overlapping meshes for the solution of partial differential 

equations J. Comput. Phys. 90 1-64 
[2] Steger J L, Dougherty F C and Benek J A 1983 A Chimera grid scheme Proc. Advances in Grid Generation 

ed K N Ghia and U Ghia (Houston) ASME FED vol 5 pp 59-69 
[3] Steger J L 1982 On application of body conforming curvilinear grids for finite difference solution of external 

flow Numerical Grid Generation ed J F Thomposon (Amsteram: North-Holland) pp 295-316 
[4] Phillips N A 1957 A map projection system suitable for large-scale numerical weather prediction J. Meteor. 

Soc. Japan 75 262-7 

[5] Phillips N A 1959 Numerical integration of the primitive equations on the hemisphere Month. Weather Rev. 
87 333-45 

[6] Browning G L, Hack J J and Swarztrauber P N 1989 A comparison of three numerical methods for solving 

differential equations on the sphere Month. Weath. Rev. 117 1058-75 
[7] Ronchi C, Iacono R and Paolucci P S 1996 The "cubed sphere": A new method for the solution of partial 

differential equations in spherical geometry J. Comput. Phys. 124 93-114 
[8] Hernlund J W and Tackley P J 2003 Three-dimensional spherical shell convection at infinite Prandtl 

number using the 'cubed sphere' method Proc. Second MIT Conference on Computational Fluid and 

Solid Mechanics (Cambridge) 
[9] Dudhia J and Bresch J F 2002 A global version of the PSU-NCAR mesoscale model Month. Weather Rev. 

130 2989-3007 

[10] Yoshida M and Kageyama A 2004 Application of the Yin- Yang grid to a thermal convection of a Boussinesq 

fluid with infinite Prandtl number in a three-dimensional spherical shell Geophys. Res. Lett. 31 L12609 
[11] McKenzie D P, Roberts J M and Weiss N O 1974 Convection in the earth's mantle: Towards a numerical 

simulation. J. Fluid Mech. 62 465-538 
[12] Turcotte D L and Schubert G 2002 Geodynamics (Cambridge: Cambridge Univ. Press) 
[13] Bercovici B, Schubert G, Glatzmaier G A and Zebib A 1989 Three dimensional thermal convection in a 

spherical shell J. Fluid Mech. 206 75-104 
[14] Kono M and Roberts P H 2002 Recent geodynamo simulations and observations of the geomagnetic field 

Rev. Geophys. 40 1013 

[15] Kageyama A, Kameyama M, Fujihara S, Yoshida M, Hyodo M and Tsuda Y 2004 A 15.2 Tflops simulation 
of geodynamo on the earth simulator Proc. ACM/IEEE Supercomputing Conference SC2004 (Pitssburgh) 



[16] Wesseling P 2004 An Introduction to Multigrid Methods (Philadelphia: R. T. Edwards Inc., John Wiley & 

Sons Ltd., 1992. Corrected Reprint) 
[17] Takahashi K et al 2004 Development of nonhydrostatic coupled ocean-atmosphere simulation code on the 

earth simulator Proc. 7th International Conference on High Performance Computing and Grid in Asia 

Pacific Region (Omiya) pp 487-94 
[18] Komine K, Takahashi K and Watanabe K 2004 Development of a global non- hydrostatic simulation code 

using yin-yang grid system Proc. 2004 workshop on the solution of partial differential equations on the 

sphere (Yokohama) pp 67-9 

[19] Takahashi K 2004 Development of nonhydrostatic coupled ocean-atmosphere simulation code Annual Report 
of the Earth Simulator for Fiscal year 2003 pp 63-7 



